####################################################################################
###
###   Immigration, Public Housing and Support for the French National Front
###   Gloria Gennaro
###
###   Paper Figure 4
###
####################################################################################


library(ggplot2)
library(multiwayvcov)

################################################################################
# Set Up
################################################################################

# Load Working Directories
source(paste0(wd_main, '/2_code/00_working_directories.R'))

# Load data
df = read.csv(paste0(wd_data, '/dataset_newspaper_analysis.csv'))

################################################################################
# Auxiliary Functions
################################################################################

LinearCombTest = function (lmObject, vars, .vcov = NULL) {
  if (is.null(.vcov)) .vcov = vcov(lmObject)
  beta = coef(lmObject)
  sumvars = sum(beta[vars])
  se = sum(.vcov[vars, vars]) ^ 0.5
  tscore = sumvars / se
  pvalue = 2 * pt(abs(tscore), lmObject$df.residual, lower.tail = FALSE)
  matrix(c(sumvars, se, tscore, pvalue), nrow = 1L,
         dimnames = list(paste0(paste0(vars, collapse = " + "), " = 0"),
                         c("Estimate", "Std. Error", "t value", "Pr(>|t|)")))
}

lincomCoeff = function(model){
  out = data.frame()
  vcv = cluster.vcov(model, cluster = df$dep)
  out = rbind(out, c(LinearCombTest(model, c("x"), vcv), nobs(model)))
  out = rbind(out, c(LinearCombTest(model, c("x", "x:factor(imm_quant_99)2"), vcv), nobs(model)))
  out = rbind(out, c(LinearCombTest(model, c("x", "x:factor(imm_quant_99)3"), vcv), nobs(model)))
  colnames(out) =  c('coef', 'se', 'tval', 'pval', 'nobs')
  return(out)
}

marignaleffects = function(model){
  out = data.frame()
  vcv = cluster.vcov(model, cluster = df$dep)
  out = rbind(out, c(LinearCombTest(model, c("x"), vcv), nobs(model)))
  out = rbind(out, c(LinearCombTest(model, c("x:factor(imm_quant_99)2"), vcv), nobs(model)))
  out = rbind(out, c(LinearCombTest(model, c("x:factor(imm_quant_99)3"), vcv), nobs(model)))
  colnames(out) =  c('coef', 'se', 'tval', 'pval', 'nobs')
  return(out)
  
}



################################################################################
# Estimation and plot functions
################################################################################

risultati = function(name){
  
  pars = as.list(match.call()[-1])
  
  # Select variables
  df$y = scale(df[,as.character(pars$name)]/df$counts_newspapers)
  df$x =  scale(df$sruT)
  
  # Main analysis
  container_lincom = data.frame()
  container_marg = data.frame()
  
  formula = 'y ~ x*factor(imm_quant_99) + factor(dep) + factor(year)'
  controls = c('', '+ factor(reg)*year', '+ factor(reg)*year + log(pop_dep)',
               '+ factor(reg)*year + log(pop_dep) + I(num_res_own/num_log) + I(num_empty_log/num_log)')

  for (ctr in controls){
    model = lm(as.formula(paste0(formula, ctr)), data=df)
    container_lincom = rbind(container_lincom, lincomCoeff(model))
    container_marg = rbind(container_marg, marignaleffects(model))
  }
  
  # Linear Combinations 
  x = scale(df$sru_weight)
  model = lm(y ~ x*factor(imm_quant_99) + factor(dep) + factor(year) + factor(reg)*year, data=df) # %>% summary()
  container_lincom = rbind(container_lincom, lincomCoeff(model))
  container_marg = rbind(container_marg, marignaleffects(model))
  
  # Clean container
  clean_container = function(container){
  colnames(container) = c('coef', 'se', 'tval', 'pval', 'nobs')
  container$group = rep(c('Low\nImmigration', 'Medium\nImmigration', 'High\nImmigration'), 5)
  container$estimation = rep(c('No Controls', '+ Regional Trends', '+ Population', '+ Housing market', '+ Weighted'),each=3)
  container$group.f = factor(container$group, levels = rev(unique(container$group)))
  container$estimation.f = factor(container$estimation, levels = unique(container$estimation))
  container_stored = container
  }
  
  
  container_lincom = clean_container(container_lincom)
  container_marg = clean_container(container_marg)
  output = list('lincom' = container_lincom, 'margins' = container_marg)
  return(output)

}


myplot = function(container, title){
  container$label = paste0(as.character(round(container$coef,3)), '\n (',  as.character(round(container$se,3)), ')')
  container$plotnum = ifelse(container$sent=='Combined', 'All Articles', 
                             ifelse(container$sent=='Negative', 'Articles with\nNegative Sentiment',
                                    'Articles with\nPositive Sentiment'))
    
  ggplot(container, aes(x=coef, y=group.f, group=rev(sent))) +
  facet_grid(.~plotnum, scales = "free") + 
  geom_point(aes(shape=sent), size=4,  position=position_dodge(width=.2)) +
  geom_errorbar(aes(xmax = coef+se*qnorm(.05,lower.tail=FALSE), xmin = coef-se*qnorm(.05,lower.tail=FALSE)), width =.2, size=.6,  position=position_dodge(width=.2))  +
  theme_bw() +
  xlab('Co-Occurrence of Housing and Immigration (SD)') +
  ylab('') +
  ggtitle("") +
  geom_text(data=subset(container, sent=='Positive'), aes(x = coef, y = group.f, label = label), size=4, vjust = 0, nudge_y = -0.4) +
  geom_text(data=subset(container, sent!='Positive'), aes(x = coef, y = group.f, label = label), size=4, vjust = 0, nudge_y = .15) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2) +

  theme(text = element_text(size=16),
        panel.grid.major.x = element_blank(),
        axis.text.y = element_text(size=12, hjust = 0.5),
        legend.title = element_blank(),
        legend.position = 'bottom') 
}


################################################################################
# Estimates
################################################################################

out1 = risultati(score)
out1_neg = risultati(score_neg)
out1_pos = risultati(score_pos)


################################################################################
# Figure 4
################################################################################

container_all = rbind(
  cbind(out1$lincom[out1$lincom$estimation=='+ Regional Trends',], sent=rep('Combined', 3)),
  cbind(out1_neg$lincom[out1_neg$lincom$estimation=='+ Regional Trends',],sent=rep('Negative', 3)))

myplot(container_all, 'Housing & Immigrants')
ggsave(paste0(wd_res, '/figures/fig4.pdf'), height=7, width=12)

